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Abstract 

The motion of finite-size particles (inertial particles) is described by the Maxey-Riley equation, which in 
its full form contains the history force. This force is represented by an integral whose accurate numerical 
evaluation is rather difficult. Here, a systematic way is presented to derive numerical integration schemes 
of arbitrary order for the advection of inertial particles with the history force. This involves the numerical 
evaluation of integrals with singular, but integrable, integrands. Explicit specifications of first, second and 
third order schemes are given and the accuracy and order of the schemes are verified using known analytical 
solutions. 
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The advection of finite-size particles (often called inertial particles) plays an important role in many 
environment-related phenomena ranging from meteorology to oceanography, e.g. cloud microphysics [1], 
as well as in engineering |2|. Particle-based modeling has been applied to the formation of planetesimals 



9< 



in the early solar system [3| and the aggregation and fragmentation processes in fluid flows Example 
applications are pollutant-transport forecasting for homeland defense [5j, and the location of a toxin or 
biological pathogen spill (e.g. anthrax) from outbreaks in a street canyon [||. Other recent results indicate 
that inertial particles might play a role in hurricane dynamics 0] and in the feeding dynamics of certain 
marine animals [8]. 

The basic equation of motion for a small spherical particle of radius a and mass m p in a viscous fluid is 
given by the Maxey-Riley equation [jjl Hoi]: 

dv Dm m f ( dv Bu\ 2 — f* 1 ( dv du\ 

m v —— = rut— — — — — bnaptv (v — u) — ba QrJirv / , — dr. (1 

p dt 1 Bt 2 \dt Bt J ^ v ' y/v Jtoy/t=T\dT dr J w 

Here, v = dr/dt is the particle velocity, u(r,t) the fluid velocity, rrif the mass of the fluid excluded by the 
particle, v the kinematic viscosity of the fluid and qj the density of the fluid. The two appearing derivatives 

du du Bu du 

—— = — — h v ■ V11 and — — = — — h u ■ vu 

dt dt Bt dt 

denote the full derivative along the trajectory of the particle and of the corresponding fluid element, re- 
spectively. The terms on the right-hand side of ^ are: the force exerted by the fluid on a fluid element 
at the location of the particle, the added mass term describing the impulsive pressure response of the fluid, 
the Stokes drag, and the history force. In this form of the equation gravity and the so-called Faxen correc- 
tions are not included. The history force accounts for the viscous diffusion of vorticity from the surface of 
the particle along the trajectory [9( and renders the advection equation to be an integro-differential equa- 
tion whose solution is much more demanding than that of an ordinary differential equation. Because of 
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this difficulty, this integral term is neglected in nearly all the applications mentioned above. Experimen- 
tal and analytic efforts [II, 12 1 indicate, however, that the history force might have significant effects for 
non-neutrally-bouyant particles in simple flows. In a recent study |13| it has also been shown that it plays 
an important role in chaotic advection of inertial particles, i.e. advection in a simple flow but with chaotic 
trajectories. The present paper will detail the derivation and analysis of the numerical schemes developed 
for the investigations in this latter study. 

The history force poses the main difficulty for a numerical integration of ((T|). There are basically three 
problems: (i) the singularity of the kernel l/y/t — r, (ii) the fact that Jl} is an implicit integro-differential 
equation due to the appearance of dv/dt on the right hand side and (iii) the high computational costs for 
a numerical integration. The first point (i) is the most involved one and will be addressed by a special 
quadratur43 scheme. The implicitness (ii) is not a major issue and can be addressed rather easily as we 
will see. The last point (iii) stems from the necessity to recompute the history force, which is an integral 
over all previous time-steps, for every new time-step. Therefore the computational costs grow with the 
square of the the number of time-steps and can become quite substantial for long integration periods. This 
difficulty is inherent to the dynamics governed by the history force and cannot be addressed without further 
approximations. Note however that a higher order scheme (typically) reduces the number of necessary time- 
steps and therefore diminishes the problem of high computational costs indirectly. Furthermore the final 
form of the numerical scheme will be formulated as a weighed sum, which is well suited for a numerical 
evaluation on modern CPU architectures. 

The correct numerical treatment of the full Maxey-Riley equation and in particular of the history force 
has received little interest in the past, in spite of an increasing number of studies supporting its importance. 
Michaelides [3] transformed the Maxey-Riley equation to a second order equation in which the history 
integral contains only the fluid velocity, but not the particle velocity. This makes the evolution equation 
explicit. Furthermore, according to Michaelides, this form of the equation allows a sparser sampling of the 
particle's history, which leads to savings in computational time and computer memory. However, the history 
integral still has a similar form as in |lj and the difficulties of an accurate numerical evaluation remain. Two 
previously proposed schemes addressing the history integral have been tested by Bombardelli et al. (la| . 
They found the accuracy of the schemes to be 0{^fh) and 0(h), where h is the time-step. In a recent work 
Hinsberg et al. [I(| have proposed a first ordei@ scheme for the computation of the history force, i.e. the error 
is 0(h 2 ), which represents a significant advancement over previously known schemes. Furthermore Hinsberg 
et al. developed a method to decrease the needed amount of history for the computation of the history force, 
by approximating the tail of the history kernel with exponential functions. This leads to significant savings 
of computational time and computer memory. This method can be viewed as a major improvement over the 
method of a window kernel where the kernel is set to zero for time lags larger then a certain window time 
[Iii . The present paper will describe the construction of arbitrary high order methods for the integration of 
particle trajectories with the history force and will give explicit specification of the first, second and third 
order methods with an accuracy of 0(h 2 ), 0(h 3 ) and 0(h 4 ), respectively. Approximate forms of the history 
kernel will not be considered. However, the developed schemes can easily be modified for a window kernel 
or the more advanced approach proposed by Hinsberg et al. 

The rest of the paper is structured as follows: First some general notes about the history force and 
the Maxey-Riley equation will given. Afterwards a numerical quadrature scheme for the history force and 
its derivation will be presented. In the next section this quadrature scheme will be incorporated into an 
integration scheme for the numerical solution of the full Maxey-Riley equation. The full integration scheme 
will then be tested against known analytical solutions. This is followed by a section on the stability properties 
of the algorithm, and by the conclusions. 



1 In this article the term "quadrature scheme" refers to a numerical scheme for the approximation of an integral whereas the 
term "integration scheme" refers to a scheme for the approximation of the solution of the whole integro-differential equation. 

2 In the paper by Hinsberg et al. the scheme is said to be of second order. This is due to a different definition of the meaning 
of "order". Here, a scheme with an error term proportional to the square of the time-step is considered to be of first order as it 
is accurate up to the first order; in the same sense as the Euler-method is a first order scheme. 
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1. Introductory Notes 



Measuring time and velocity in units of T and U, the dimensionless Maxey-Riley equation becomes 



1 dv Dm 1 , , 3 1 r 1 /du dw 



^ - -jr ( 2 ) 



i? dt Dt Si v ' V 7T St J to sJt^T \ dr dr 
Here two dimensionless parameters appear, the density ratijl 

mf + zrrip 

and the Stokes number 

st- 1 -^ 

bt ~3 T ' 

a ratio of the particle's viscous relaxation time and the characteristic time of the flow T. In smooth large- 
scale flows there is often only one typical time scale whereas in a turbulent flow there a many. In the latter 
case the smallest time scale is appropriate and the Kolmogorov time is a frequent choice. 

An important condition for the validity of equation ([2]) is that the particle Reynolds number Re p = 
\v — u\ ajv remains small during the entire dynamics 0. In addition, the Stokes number must be small 
(i.e. the particle's typical time scale is much smaller than that of the flow) and a -C L. The last condition 
assures that the so-called Faxen corrections are negligible 

Many of the derivations and concepts in this article are applicable for any kernel appearing in the history 
force integral. Therefore, in the following, a general kernel K (t — r) will be used where the derivations do 
not depend on the particular form of the kernel. The explicit specification of the quadrature scheme and 
the tests of the numerical schemes will be given for the standard kernel 

K{t-r) = ^=. (3) 

Before we proceed with the derivation of the quadrature scheme, let us first rewrite the history force 
integral in a different form 



f K (t - r) A /(t ) dr + K(t - t )f(t ) = ± [ K(t- 



i/(r)dr, (4) 



where /(r) = v — u. This relation can be verified using integration by partfl Equation ([TJ) has been derived 
with the assumption of a particle starting with the same initial velocity, as the fluid, i.e. v(to) = u(to). In 
this case the second term on the left-hand side of ^ vanishes. In the case of different initial velocity the 



additional term (tt(io) ~ v (to)) /V^ ~ has been given in [14J and [17|, which is exactly the additional term 
appearing in ((U). Therefore the Maxey-Riley equation can be written in the following form, which is now 
also valid for initial conditions with v(t§) ^ u(ta), 



1 dv Dm 1 . . 3 

= (v — u) ~ \ — 

Rdt Dt St y ' V it St dt 



r) (y — u) . (5) 



It is beneficial to use this form of the history force (the rhs of because it enables us to compute an 
integral of the history force by dropping the derivative. This improves and simplifies the numerical scheme 



3 In some cases the density ratio is defined as R = 2m^/(mj + 2m p ), which differs by a factor of 3/2 from the definition 
here. 

4 When the kernel has singularities, one has first to use integrals with the upper bound of t — e, then perform integration by 
parts and finally take the limit e — > (to prevent the appearance of singularities outside of integrals). An alternative for the 
standard kernel is to use a transformation of the integration variable r — > x = \/t — r. 
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as we will see. At this point it is interesting to note that for the standard kernel the history force is equal 
to a fractional derivative of the Riemann-Liouville type: 

Thus the numerical methods developed here can be also considered as high order methods for the numerical 
computation of fractional derivatives and the solution of fractional differential equations. 



2. The Quadrature Scheme 

In this section a systematic way is presented for the construction of quadrature schemes for integrals of 
the type 



f K(t-T)f(T)dr. 

J t 



When the kernel K is a, well behaved function no special effort is needed and standard schemes can be used. 
However in cases where the kernel has an integrable singularity, like e.g. the standard kernel Q at r = t, 
standard numerical methods, like the Newton-Cote^l schemes, lead to large errors as we will see. This is 
due to the necessity to evaluate the whole integrand including the kernel near the singularity. We will avoid 
this by constructing a specialized scheme in which the kernel is already integrated analytically. 

Due to the linearity of the history integral with respect to / any quadrature scheme for this term can be 
expressed as a weighted sum 



#(t-T)/(r)dT«5>j/(Ti), 

-° 3=0 



where n = to + hi, n — it — to) /h and h is the time-step. The main topic of this section is the derivation 
and specification of the coefficients /ij. The general procedure is to first split the integral into intervals of 
length h 



ft n 1 rri+i 

K(t-r)f(r)dT = J2 K{t-r)f{r)Ar, 

J t t=Q J Ti 



then to approximate /(r) in every of the intervals with a polynomial and finally to compute the appearing 
integrals analytically. The order of the polynomial will determine the order of the scheme. 

Let us first examine the simplest case of a linear approximation, i.e. an order one scheme. By approxi- 
mating /(t) linearly in the interval [r.;, Tj + i] 

/(r) = f( Ti ) + /(Tm) fe ~ /(Ti) (r - n ) + O(h') (6) 

we obtain 

£ +1 K(t-r) /(r) dr = (f( n ) + 0{h 2 )) £ K (t - n - r) dr + /(T ' +l) / ~ /(Ti) jT tK (t - n - r) dr. 
In many cases the appearing integrals can be computed analytically, e.g. for the standard kernel 



rn+i 



/to t*,-s , r o u _ - iK f(n+i)-f(n) 



■ dr = {f{n) + 0(h 2 )) [-2Vt^F~ 



r 



4, 3 



-2rVt-Ti - t - - (t - n - r) 



h 



5 Well known Newton-Cotes schemes are e.g the trapezoidal rule and Simpson's rule. 

6 The error of an approximation will be denoted by 0(h m ), i.e the error is bounded by Ch m with a fixed constant C . 
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Summing up the terms for each of the intervals we obtain a formula for the whole integral, e.g. for the 
standard kernel 



=<lr = 2f(t )Vt=h + f{n+1 \ f{T%) { {t ~ T ' +l)l ) +0(h 2 )Vt^h- (7) 

i=0 



/to V* 



Note that no singular or diverging expressions appear. For this it is crucial to approximate only f(r) with 
polynomials, but not the whole integrand. 

As already mentioned the quadrature scheme is linear in / and can thus be expressed as a weighted 
sum. Such a form is best suited for a numerical evaluation as modern processors/compilers can optimize 
this kind of operations rather well. We will index the coefficients of the sum in reversed order, i.e. we use 
the sum . Hjf(r n -j) instead of y\ fJ*jf(Tj). This is more natural as it turns out that the coefficient of 
/(tj) depends on n — j. For the standard kernel the coefficients for the first order quadrature scheme can 
be obtained from 



dr = VhY,aqf(T n -j) + 0(h 2 )Vt^h 




3=0 

| <! (i - 1) 3/2 + (i + 1) 3/2 - 2i 3 / 2 0<j<n (8) 



(n - if 2 - n 3 / 2 + 



3 = n. 



Here the factor y/h has been pulled out of the coefficients to make them independent of the time-step h. Also, 
note that the coefficients aj depend on n, the number of intervals for the approximation of the integral. The 
first order scheme specified by (J7J and © is equivalent to the one presented in , although the equivalence 
is not obvious. 

The procedure just shown can be generalized to derive quadrature schemes of arbitrary high order. The 
basic ideas stay the same, however the technical details make the derivation complicated. Here, only a 
simplified overview of the construction will be given. The full derivation with all the technical details can 
be found in |Appendix A| 

To obtain a quadrature scheme of order to, we approximate / in every interval [Tj,Tj_|_i] with an m-th 
order polynomial and solve the remaining integrals analytically. An interpolating polynomial of order to 
is uniquely determined by the values of / at m + 1 time-points. Let us denote these time-points by 0^, 
where i is the index of the interval and k 6 {0, . . . , to}. Using the Lagrangian representation of polynomial 
interpolation we obtain the approximation in the i-th interval 

TCI TCI * 

f(r) = £ f(e ik )L ik (r) + 0(h m +') L lk (r) = J] f^J-. 
fe=o i=o ° ik Ua 

The time points 6^ can in principle be chosen arbitrary. However, it is clear that this choice will strongly 
influence the quality of the interpolation. Obviously, the points Tj and r; + i should be included when 
interpolating in [Tj,Tj+i]. These time-points were our choice for the first order approximation ©. For 
higher order approximations we need more points additionally to r% and Tj_|_i. Reasonable choices are the 
points closest to the bounds of the interval, i.e. Tj_i, Ti + 2, ^-2, ■ • • (given we want to stay on the grid defined 
by the Tj). And indeed we will use t.;_i, Tj, Tj + i for the second order approximation and Tj_i, Tj, Tj + i, r^ + 2 
for the third order approximation. This can be generalized to arbitrary orders by choosing Oik = T i—\ m/2\+k 
where the operation |_-J denotes taking the integer part, often called the floor function. 
With this definitions we can express the history integral as 



„t n— 1 m rT+i n—X m 

/ dTK(t-T)f(r) = 1 £Y,f( e ^ drK(t-r)L k (r)+E = Y / T,- 

to i=0 k=Q izi . i=0 k=0 



where E = f. K(t — t)<1t 0{h m+1 ) is the error term. We naturally obtain a weighted double-sum due to the 
use of the Lagrangian representation of the interpolating polynomial, where the f(9ik) appear as coefficients 
of the polynomials Lj/.. Compare this with derivation of the first order scheme where we started from the 
linear interpolation ([5]), which is not in the Lagrangian form and thus a reordering of terms was necessary 
to get from © to ©. 

The integrals A^fc do not involve /(r) and can be computed in advance; for many kernels even analytically, 
including the standard kernel. Now the final step is to reorder the double sum to a single weighted sum 



L 



I n — 1 rn n 

dr K (t - r) f(r) = £ £ f(e ik )X ik = £ fiffa-j). 

*° i=0 k=0 j=0 



This procedure is detailed in |Appendix A| Note that the coefficients /x" (just like a") have a dependence on 
n (see |Appendix A| for details). 

In the following the coefficients for the standard kernel ^ will be given for a second and third order 
scheme, denoted by /3™ and 7™ respectively. The factor \fh has been extracted from the coefficients so that 
they do not depend on the time-step h. 

The second order approximation is 

Jtn Vt — T ._„ 



with /3™ for n = 2 and n = 3 



A/=0,l,2,3 



- ^ 



and for n > 4 



|V2 



A ( (j + 2 f 2 - 3 + 1) 5/2 + 3// 2 - (j - I) 5 / 2 ) 
+ | (- (3 + 2) 3/2 + 3 (i + 1) 3/2 - 3// 2 + (j - I) 3 / 2 ) 

^ (-2n 5 / 2 + 3 (n - 1) 5/2 - (n - 2) 5/2 ) 
+ §(4n 3 / 2 -3(n-l) 3/2 + (n-2) 3 / 2 ) 



^ (n 5 / 2 - (n - 1) 5/2 ) + § (-3n 3 / 2 + (n - f ) 3/2 ) + 2^ j 
The third order approximation is 

rt 1 n 

/ dr ^=/(r) = VhJ2 7" /(r„-i) + 0(fr 4 )V*^ 



3=0 
j = 1 
J = 2 

2 < j < n - 1 



j = n - 1 



fi 



with 7™ for 3 < n < 6 



7j 5 =o..5 



7j=o.. 



3 68 r- 6 r- 12 R 16 /- 

7j=o..3 = T^V3; ^3; -V3; — V3 

244^ 1888 _ 976 ^ _656 488 ^ 544 _ 976 ^ _292 244 ^ 

315 ' 315 315 ' 105 105 ' 105 315 ' 315 315 

V2; ^73-^^2; ^VS-^Vs+^Vl; ™V3-™V2; 

'105 315 '63 105 105 ' 21 35 315 ' 



4 _ 244 



244 
315 
220 



310 1U5 315 

220 r- 1448 /- 244 /- 
IT^" 105-^+315^ 



and for n > 7 



244 

315 
1188 
35 



164 
"63 



V5 



, 362 /- 
+ 105^ 



1UO OlO OO 1UO 

V2; ^-^;^-^^+^V2; ^ - ™ + ™ 
105 315 315 105 105 21 315 35 



J ' 105 315 315 105 
^ 11168 _ 1448 ^ 244 
105 105 



105 105 21 315 35 3 

244 /- 936 k 22336 362 - 754 r- 5584 

H V2; V6 1 V3; VfH 

315 35 315 105 105 315 



y3-^V2; 
315 ' 

558 
+ 315 



362 /77 976 /TT. 

105 : l ■ x 

5584 

315 105 v " 1 105 



1* = \ 



315 

ii«y3+|§|V2 



1130 /c _ 22336 , 724 fn 976 fn 
~63~ V ° ~315~ + ~35 V 315 VZ 



^ ( (j + 2) 7 / 2 + (j 2) 7 / 2 - 4 (j + 1) 7/2 - 4 (j - I) 7 / 2 + 6/^ 



+ f 



4 (j + 1) 3/2 + 4 (j - 1) 3/2 - (j + 2) 3/2 - (j - 2) 



3/2 



n3/2 



\3/2 



•3/2 ' 



(n 7 / 2 - 4(n - 2) 7/2 + 6 (n - 3) 7/2 - 4(n - 4) 7/2 + (n - 5) 7/2 ) - ^7 



+|- 3/2 



+ § (n - 2) 3/2 - | (n - 3) 3/2 + | (n - 4) 3/2 - f (n - 5) 3/2 



((n - 4) 7 / 2 - 4(n - 3) 7/2 + 6 (n - 2) 7 / 2 - 3n 7 / 2 ) + §n 5 / 2 
-2n 3 / 2 - | (n - 2) 3/2 + § (n - 3) 3/2 - § (n - 4) 3/2 




2) 



7/2 



+ (n - 3) 7/2 ) - fn 5 / 2 + 4n 3 ^ 2 + f (n - 2) 3 / 2 - f (n - 3) 3/2 
(n - 2) 3/2 + 2^n 



22 3/2 



3 = 

i = i 
i = 2 
i = 3 

3 < j < n - 3 



j = n - 



n - 
n. 



Now that the quadrature schemes are fully specified, let us verify the correctness of the derivation and 
in particular the order of the schemes by using a test case where the analytical solution of the integral is 
known. We choose the case /(r) = sin(r) where the history integral can be computed with the help of the 
Anger function J v {t) [18], which is a generalization of the Bessel function J n [t) to fractional values of n, 



^ J==LdT = (j, (t) - J_l (*)) = I{t). 



(9) 



To verify the order of the scheme let us analyze the global error 

e(h) = max \I(t) - I num (t, h)\ , 
te[o,io] 



where I(t) denotes the exact value of integral in Q and I num (t, h) the numerically approximated value. 
Fig. [1] shows the dependence of this global error on h for the three numerical quadrature schemes given here 
(specified by a™, /?", 7™) and a second order, semi-open Newton- Cotes scheme jTsj- We see that errors of 
the schemes are proportional to h m+1 for the m-th order scheme, thus verifying the order of the quadrature 
schemes (at least for this test case). Also we see that a standard second order quadrature scheme (the 




Figure 1: Scaling of the global error e(h) of the quadrature schemes for the test-case /(r) = sin(r). 



Newton-Cotes scheme) performs very badly as the error scales only with yh. This is also true for higher 
order Newton-Cotes schemes and is due to the necessity of a numerical evaluation of the kernel near the 
singularity. 

The correctness of the quadrature schemes has also been tested using the analytically treatable case of 
a polynomial of arbitrary order and led to similar results. 



3. Integration of the full Maxey-Riley equation 

In this section the quadrature scheme developed in the previous section will be incorporated in a multi- 
step integration scheme for the full Maxey-Riley equation. To this end we formulate the Maxey-Riley 
equation for the velocity difference w = v — u in a given flow held u: 



dw ,„ , du „ „ R „ /3d /"* . , x , 

- = (R - 1) - - R W • Vu - - W - R^ — - j K(t - r) W (r) dr. 



(10) 



Together with the evolution equation for the particle position 

dr 

— = v = w + u 
dt 

equation (fT0|) fully specifies the motion of an inertial particle in a fluid. Integrating (fTU)) from t to t + h and 
using the abbreviations 

G = (R - 1) ^ - Rw ■ Vu - ^-w 

dt St 



we obtain 

(•t+h 



G(T)dr + H(t + h)- H(t). (11) 



Here the integration of the history term can be performed trivially due to relation This simplifies 
the integration scheme considerably. Furthermore, we now have to compute a history integral of w and 
not dw/dr, where the former will generally fluctuate slower and is therefore better suited for a numerical 
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quadrature. The history integral H can be computed with the schemes developed in section [5] The integral 
of G can be approximated using polynomial interpolation. We use only the present and previous values of 
G for this approximation in order to obtain an explicit scheme: 



/i-t-n 
G(r) dr = hG(t)+0(h 2 ) 

rt+h u 

/ G(r)dr = -(3G(t)-G(t-h))+0(h 3 



rt+h h 

/ G(t) dr = — (23G(i) - 16G(i - h) + 5G(t - 2h)) + 0{h A ). 
Jt 12 

These expressions can be found in any literature on Adams-Bashforth multi-step methods, e.g. [l9| . 
A final point which we have to consider before writing down the complete scheme, is that 

n+l 

H(t + h) = Yl ^ +1 w(r n+1 - J ) + 0(h m ) 

3=0 

depends on w(r n+ i) = w(t + h) and thus can not be evaluated before w(t + h) is known. This is due to the 
implicitness of the Maxey-Riley equation. However this is easily dealt with by bringing w(t + h) to the left 
hand side of (fTTI) . 

If we now consider (fTTj) on the grid t n = <o + nh, define £ = Ry/3/(nSt)^/h and use abbreviations of the 
type w n — w(t n ) we can specify the complete integration schemes of first, second and third order for the 
Maxey-Riley equation: 



n+l 



h{w n + u n ) +0{h 2 ), 



(1 + W n+1 = W n + hG n -^J2( a l+l W n-3 - oQWn-j) + VU~hO(h 2 ), (12) 

3=0 

r n +i =r n + - (3 (w n + u„) - (tu„_i + tt„_i)) + G(h 3 ), 

h n 

(1 + ^o +1 ) = w n + - (3G„ - G„_r) - ^ ^ (jfi^i to n -j - + VV^O^ 3 ), (13) 

r„+i = r„ + ^ (23 (io n + tt„) - 16 (u>„_i + u„_i) + 5 (tu„_2 + Wn-2)) + 0{h A ), 

h n 

(1 + £ 7 ? +1 ) = ^ + 12 (23G ™ " 16G «- X + 5G "- 2) - ^ £ (77+>n-, - T>„-i) + ^/t~r a o{h A ). 

3=0 

(14) 

The coefficients a", /3™ and 7™ are given in section [2] One reason to include the first and second order 
schemes here (besides the third order one) is that one cannot start the integration with the third order 
scheme as the previous values G n _i, G n _j are not available at the beginning. This is a problem common 
to all multi-step methods. The simplest solution is to use the first and second order schemes for the first 
two steps and the third order scheme for the rest. To perform the first step of the integration (n — 0) with 
(fT2"]l the coefficients a° are needed, which we define to be zero as no history is present at t = to- Ideally, 
we would perform the second step (n = 1) with (1131) . But /3" is defined only for n > 2, leaving us with 
two options: (i) perform the second step with the first order scheme or (ii) define m = a 1 and accept a 
reduced accuracy. The second option is at least as accurate as the first one and will thus be assumes in the 
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Figure 2: Procedure to start the integration with multi-step methods. 
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Figure 3: (a) The exact trajectory of a particle starting at ro = (1, 0) with too = and the parameters R = 0.75 and St = 0.3. 
The dots show the position at integer time units. Also shown are the approximations of first, second and third order for 
h = 10~ 2 , where the latter two are overlapped by the exact trajectory and are thus not visible, (b) The relative error of the 
the numerical solutions obtained by the first, second and third order schemes IU2I - II14H with h = 10~ 2 . 



following. The same considerations are applicable to the third order scheme ([HI) , leading to the definition 
7? = /3? and allowing its use for n > 2 (instead of for n > 3). 

Using lower order schemes for the first two steps makes them less accurate. A more advanced strategy 
is to begin the integration with a smaller time-step to account for the reduced accuracy and switch on the 
third order scheme with the normal time-step when it is applicable. This procedure is demonstrated in 
figure [21 At the beginning of the integration, eight small steps with time-step h! = h/A are taken. From 
the time-point to + 2h on the third order scheme can be applied with the normal step size h as enough 
previous values are available then. In practice the size h! of the small steps can be much smaller than h, 
e.g. h! = h/100. This procedure can be further refined, e.g. in figure [5] it would be sufficient to take steps 
of h/2 in the interval [to + h, to + 2h]. However, the savings in computational time due to this optimization 
will generally be not large enough to compensate for the increased complexity of the algorithm. 

The integration methods (p~2|l - (fl4)) can be viewed as an extension of the Adams-Bashforth multi-step 
methods to the case of an integro-differential equation with memory. In its present form the quadrature 
schemes in section [5] are best suited for multi-step methods with a fixed time-step. In other schemes, for 
example Runge-Kutta, half-steps are necessary, but they cannot be evaluated with the current formulation 
of the quadrature schemes. Furthermore, multi-step methods allow to profit from the fact that an integral 
of the history force can be computed by simply dropping a derivative (see (llip and comments below). 

To test the accuracy of the whole integration scheme, the motion of a particle in the flow u(r) = \r\ e v 
(rigid body rotation) will be considered. Fortunately, in this case there is an analytical solution for the 
full Maxey-Riley equation found by Candelier et al. [12J. Qualitatively, the solution is a spiraling motion 
outwards or inwards depending on whether the density of the particle is larger or smaller then that of the 
fluid, i.e. R < 1 or R > 1. Asymptotically the distance of the particle from the center grows exponentially, 
\r(t) \ ~ exp(Ai). The ejection rate A depends on the presence of the history force and thus the trajectories 
of particles with memory and without memory deviate rather quickly. This makes this flow a good choice 
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Figure 4: Scaling of the global error e(h) as a function of the time-step h. The particle was started with the initial condition 
ro = (1, 0) and wq = and the parameters R = 0.75 and St = 0.3. 



for a test of the integration scheme as an inaccurate computation of the history force is expected to lead to 
strong deviations from the analytically known trajectories. 

Figure [3^i shows the exact solution together with the numerical solutions of first, second and third order 
obtained by p^|) - (TH)l with h = 1CP 2 . Only the the first order approximation is visible, whereas the second 
and third order ones are overlapped by the exact trajectory. To understand this let us examine the relative 
error 

P /, t \ \r(t) -r num (t,h)\ 
E ^ h ) = ' 

where r num (t,h) is the numerical and r(t) the exact solution. Figure \3jp shows this quantity for h — 10 -2 . 
We see that the error improves by approximately two orders of magnitude for each additional order of the 
scheme, thus explaining the overlapping of the second and third order approximations by the exact solution 
in figure [3^,. Figure \3jp also gives information about the quality of the approximation as a function of time. 
For example, at t = 100 the first order approximation has a very large error of ca. 60%, whereas the second 
and third order approximations are rather accurate with errors of ca. 0.4% and 0.003%. At t = 100 the 
distance of the particle from the center is |r(100)| ~ 31 whereas for a particle without memory (i.e. when the 
history force is neglected) it is «476, showing that the history force has a strong influence on the particle's 
motion. Thus an accurate computation of the history force is essential for a high precision approximation, 
as obtained by the second and third order schemes. 

To examine the dependence of the error on the width of the time-step let us again use the global error 

e(h) = max \r(t) - r num (t, h)\ , 
te [o,ioo] 

where r num (t, h) is the numerical and r(t) the exact solution. Figure 0] shows that the error of the m-th order 
scheme scales as h m . This is consistent with the one-step error 0(h m+1 ) of the numerical schemes (fTS]) - (fT¥j) 
as the global error is also proportional to the number of time-steps. Thus the global error is expected to 
behave as t max /hO(h m+1 ) = 0(h m ), where t max is the integration length. 

From the dependence of e on h we can see that rather small global errors are achievable with moderately 
small time-steps when the second or third order scheme is used. 

4. Stability of the Integration Scheme 

An important property of numerical algorithms is stability, i.e. errors remain bounded during the 
iteration of the algorithm. For ordinary differential equations numerical stability is usually determined by 
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order of the scheme 


1 


2 


3 


with memory 


4.7627 


0.9428 


0.3886 


without memory 


2.0000 


1.0000 


0.5455 



Table 1: Stability thresholds h t h of the numerical schemes IU2I - H14H compared with those without the history force in HIGH . 



applying the integration scheme to the equation 

t = ^ (15) 

and verifying whether the numerical solution converges to zero. To check the stability of our schemes we 
use the equation 

dr = -H*° + s/ tI ^ dT J' ,16) 

which is the Maxey-Riley equation ([5]) in still fluid (it = 0) with R — 7rfc/3 and St = ir/3. The solution of 
this equation converges to zero algebraically (~ i -15 ) in contrast to an exponential convergence for (fT5j) . 
In general A: is a complex number. However, here the analysis is restricted to purely real and positive values 
of k. In this case, k can be set to 1 by rescaling the time and we can analyze stability as a function of the 
time-step h only. 

For the case with memory an analytical investigation of the stability is rather difficult because the linear 
recurrence relation, which results when applying the numerical scheme to the integro-differential equation 
(TT6"]) . contains every previous time-point. To determine the stability, one needs to find the eigenvalues of 
this recurrence relation. It is not clear whether this is tractable at all with purely analytical methods. Thus, 
the investigations here are restricted to an "experimental" stability analysis: the scheme is iterated for at 
least 10 6 time-steps and it is checked whether w n converges to zero. This procedure has been carried out 
for a large number of different values of the time-step h and it has been found that w n either converged to 
zero or became infinite. These two regimes are separated by the stability threshold h t h, i.e. for h < /ith 
the iterated scheme converges to zero and is thus stable and for h > hth the iterated scheme diverges and 
is thus unstable. Table [TJ shows the stability thresholds (determinde by the bisection method) for the first, 
second and third order schemes (|12l) - (fl"4"]) . For comparison the row "without memory" contains the stability 
thresholds of the schemes without the history force, i.e. normal Adams-Bashforth scheme^- For the first 
order method the inclusion of the history force increases the stability threshold, i.e. the scheme becomes 
more stable. In the case of the second and third order schemes the inverse is true; the stability threshold is 
slightly lower and the schemes are less stable when memory is included. However the influence of the history 
force on the stability of the schemes seems to be rather weak as the stability thresholds vary only by a factor 
of order unity. Summing up, one can say that the integration schemes (|12[) - (|14|) seem to have very similar 
stability properties as the corresponding Adams-Bashforth methods for ordinary differential equations. 



5. Conclusion 

This paper presents a systematic way to derive higher order numerical integration schemes for the full 
Maxey-Riley equation including the history force. Due to the singularity of the integrand of the history force 
a special numerical scheme is needed. Explicit specifications of the numerical schemes of first, second and 
third order with an accuracy of 0(h 2 ), 0(h 3 ) and 0(h 4 ), respectively, have been given. Furthermore the 
correctness and the order of the schemes have been verified by comparison with known analytical solutions. 
The accuracy of the second and third order schemes represents a substantial improvement compared with 
the methods available in the literature. 



7 The computed thresholds without memory are consistent with the known stability regions of the Adams-Bashforth methods. 
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Up until now the lack of higher order integration schemes for the full Maxey-Riley equation has been 
one of the major hurdles for research concerning the history force. I hope that the schemes developed in 
this paper will resolve some of the hurdles and thus facilitate investigations on the relevance of the history 
force in the motion of inertial particles. 

Acknowledgment 

I would like to thank Tamas Tel and Michael Wilczek for useful discussions. 

Appendix A. Details on the derivation of the quadrature scheme 

In section[5]the derivation of the quadrature scheme was presented in a simplified, not fully detailed way. 
This appendix will present the technical details and give a complete, but somewhat laborious derivation. 
To interpolate /(t) in the interval [Tj,Ti_|_i], Lagrangian polynomial interpolation is used: 

rn m 

f(r) = £ fiODLT (r) + 0(h m+1 ) L^{r) = J] ■ 

k=0 1=0 



ik °il 



Here the full dependence of the time points Of™ on n and m has been written out explicitly. In section [5] we 
have chosen #™ fc m = Ti_ i m/21 +fc - Thus the dependence on m and i is obvious, and we will see in a moment 
why a dependence on n is necessary. The problem with the above definition of 6"™ is that for i < [m/2\ we 
would obtain time points outside the integration interval [to, t], e.g. (9q™ = 7"_ 1 m/2J = ta~ L TO /2J h, and thus 
would have to rely on values of f(r) that are not available. A similar problem occurs for i > n— m+ \ m/2\ . 
To solve this, we need a definition of 9™™ that deals with the special cases i < [m/2\ and i > n — m+ [m/2\ . 

is 

'0 < i < [m/2\ 

i - [m/2\ [m/2\ < i < n - m + \m/2\ 
n — m n— m + [rn /2J < i < n — 1 

and 6*™ fc m = T ™ +fc . The offset is defined such that it is equal to i — [_m/2\ (corresponding to our naive 
ansatz = T i _y m / 2 \+k) where possible and is set to and n — m where we would obtain time points 
outside the integration interval [to , t] . 

The interpolating polynomial for /(r) in the interval [ri,Tj-)_i] can now be expressed as 

m m 

X T 1 1 TT T — T nm,] 

x \ rum /_\ 1 /-/-(/;_m+l\ r nm l _\ I u ih. ^ L 



f(r) = £ f(r 0?m+k )LT (r) + 0(h m+1 ) L^{r) = \{ 

1=0 



fc= 1=0 r °?r+ fe r °rr+' 



and integrated to yield 

drK(t- r) f(r) = £ £ f(r 0?m+k ) / dr K (t - r) £"™(r) = /(r„ ? « +fc )\T + B . 

i=0 k=0 ZZl ' i=0 k =° 

(A.l) 

where E — f tg K(t — r)dr 0(h m+1 ) is the error term. 

Let us now reorder the double sum to a single sum of the type /z™ m /(r„_j). As already mentioned 
in section [21 it turns out as beneficial to index the coefficients u™ n in reversed order, i.e. /Xg m and /x™ m 
correspond to f(r n ) and /(to) respectively. For the following calculations we will use the theta function, 
which is defined here in the following way: takes logical conditions as arguments and has the value 1 if 
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the condition is satisfied and otherwise, e.g. 8(i < 0) is equal to 1 when i < 0. The double sum in (|A.1I) 
can be expressed as a single sum 

n — 1 m n n 

E E f(r ^ +k )xT = E E 9 + fc = n - A ™ = E ti m f(T»-3)> 

i=0 fc=0 j=0 i,k j=0 

with the coefficients 

n—1 m 

= E E 9 (°" m + k = n - 3) A "™- 

i=0 k=0 

Using the definition of o" m , the sum over i can be split into three terms (for the purpose of a compact 
presentation the indices n and to will be omitted and the abbreviations a — \_m/2\ and b = to — \_ni/2\ will 
be used): 

a— 1 m n—b m n — 1 m 

f^j = y^ J /^ & (k = n - j) X ik +} j }^Q (i - a + k = n - j) Xjk + E {n - m + k = n - j) \ ik . 

j=Q fe=0 i=a fc=0 i=n-b+l k=0 

The conditions in the theta functions can be used to get rid of one summation. For example in the first 
term the condition k = n — j is satisfied at most for one value of k and thus k can be replaced by n — j. 
However one has to keep in mind that the condition may be not satisfiable at all (it is satisfiable when 
< n — j < m). Applying this kind of reasoning to the other two terms yields 



/Ltj = 6 (0 < n — j < m) h. n -j + ^6(a<n-j-fc + a<n-6) \n-j-k+a,k 

i=0 k=0 

n-1 

+ e (o < m - j < m) a * 



m—3- 
i—n — b+1 

In the second term the summation over i has been removed. The satisfiability condition depends on k and 
thus has to remain inside the sum. Simplifying the conditions we obtain 

a — l m n—1 

Hj = 6 (n - to < j < n) At, n -j+y^ 6 (to - j < k < n - j) X„- 3 -k+a,k+Q (0 < j < to) ^ Xi, m - V 

i=0 k=0 i=n-fe+l 

The condition in the second term can be used to narrow the summation range, and we obtain the final 
expression for the coefficients jj^ m 



a— l min(m,n-j) n—1 

= e(n- m <j<n)J2K:n- 3 + E K-i-k+a.k + © (0 < j < m) ^ X?™_j. 

i=0 fc=max(0,m— j) i=n — b+l 

For the case of the standard kernel ((3]) the integrals A™ can be computed analytically and thus the 
coefficients /i" m . This has been done by means of the computer algebra system Maple and the resulting 
expressions for the coefficients are given in section [2] 
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